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1 Abstract 

An overview of the current status of time depen- 
dent algorithms is presented. Special attention is 
given to algorithms used to predict fluid actuator 
flows, as well as other active and passive flow' control 
devices. Capabilities for the next decade are pre- 
dicted, and principal impediments to the progress of 
time-dependent algorithms are identified. 

2 Introduction 

Continuously expanding computer capabilities al- 
low more attention to be devoted to the simulation 
of unsteady flows. At the turn of the millennium, 
practitioners routinely compute complex 3-D steady 
flows involving 10 b — 10' grid points, and ‘2-D un- 
steady flows involving 10 5 — 10 l:i points. These feats 
are performed while carrying five or more variables 
per node! If Moore’s law persists (a fixed cost dou- 
bling of computer resources every 1.5 years) the next 
decade will provide practitioners with the resources 
to routinely simulate 3-D unsteady flows on lO 0 grid 
points. This computer capability will enable the 
burgeoning field of aerodynamic flow control (AFC), 
which is often time-dependent. 

Active flow' control offers the aerospace commu- 
nity the opportunity to expand the flight envelope 
through the use of steady suction/blowing, zero net 
mass synthetic jet actuators, or pulsed jets. These 
flow control devices exhibit promising flow control 
capabilities including separation control, thrust vec- 
toring, mixing enhancement, noise control, and vir- 
tual shape change. Benefits of flow control in- 
clude reduction in part-card count, empty weight, 
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manufacturing costs, operating cost, fuel burn, and 
noise. A number of active flow control concepts have 
been tested in the laboratory and flight. Examples 
include leading-edge suction for transition delay, 93 
zero net mass separation control 123 ' 124 ’ 125 ' 126 ' 12 ‘ 
and thrust vectoring fluidic injection. 111 Computa- 
tional studies have demonstrated that Reynolds av- 
eraged Navier-Stokes (RANS) methodologies pro- 
vide qualitative insight into active flow' control ap- 
plications. However, quantitative agreement is lack- 
ing between the computational and experimental re- 
sults. To get from the bench-top to real applications 
of flow control, reliable computational fluid dynam- 
ics (CFD) design tools must be developed and val- 
idated with the experimental and flight databases. 
An extensive amount of research is still needed to 
develop a production-type tool for active flow con- 
trol applications for the design engineer. 

A critical assessment of the current capabilities of 
time-dependent CFD, and identification of impedi- 
ments that still exist is timely. We focus on identi- 
fying the critical areas (algorithmic and modeling) 
that possess notable leverage to the success of 3-D 
AFC computations. 

The review of this material will be presented with 
the following strategy. Each section will begin with 
a broad overview of current state of the art in that 
field, follow'ed by a description of general bottle- 
necks, and specific impediments for time-dependent 
AFC computations. Finally, each section will con- 
clude with a brief summary of NASA Langley Re- 
search Center’s (LaRC) present research aimed at 
alleviating the bottlenecks, recognizing that some 
impediments are not being addressed due to limited 
resources. The fields of CFD and turbulence model- 
ing are nearly boundless! To limit the scope of the 
review, only those methodologies which have show'n 
promise in AFC simulations will be addressed. 

The paper focuses on the general areas of algo- 
rithmic issues and turbulence models, and on the 
specific area of fluid actuators. The paper is orga- 
nized as follows. Section 3 describes discretizations 
in space and time. The section begins with a broad 
discussion of the advantages of high-order schemes, 
followed by specific discussions of temporal and spa- 
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tial discretizations. Section 4 describes algorithmic 
considerations related to convergence acceleration. 
Section 5 describes the current state of turbulence 
modeling for time-dependent flows. Section 6 de- 
scribes specific considerations for effective actuator 
boundary conditions. Section 7 presents conclusions. 

3 Discretizations: Time and Space 

3.1 Wh y H igli-Ordcr? 

For reasons of efficiency, high-order schemes have 
long been advocated for use in time-dependent prob- 
lems. In 1902, Kutta recognized the virtues of inte- 
grating ordinary differential equations (ODEs) with 
high-order schemes. The following simple example 
illustrates this point. Local error e; committed dur- 
ing one step of a temporal integration is described 
by the formula 

lie, -|| < (At)f 1 (1) 

where p is the temporal order of the integration for- 
mula. Global error at time Tf is estimated by sum- 
ming all local errors after transporting each to the 
final time Tf. Estimates of global error, though not 
sharp, are generally expressed in the form 55 

Ill’ll < (±t) p max C{exp[L(T f -T 0 )]-l) (2) 


Assume a uniform grid xj = j Ax with Ax = ‘2 tt/N. 
A general «; T n r + 1 point spatial discretization is 


1 " r 

U x \ j = ^ ^2 a iU(xj+i), j = 0,N-1 ( 5 ) 

l= — ni 

Substituting equation (5) into equation (3) and solv- 
ing the system of ODEs ( 'Ut + MU = 0) in Fourier 
space yields the modal solution 

U(x,t) = c < * (* - *(*) t) ( 6 ) 


where a(k) is the wave speed of the semi-discrete 
problem and is related to the Fourier image of the 
spatial operator. For example, the second- and 
fourth-order central difference waves speeds are 


«(*)2 


a{k) 4 


sin(kAx) 

£Q> - ■ ■ 

2 k Ax 

8 sin(kAx) — sin(2kAx) 
' ° 12 kAx 


( 7 ) 

( 8 ) 


For real a (/«•), the difference (error) between the ex- 
act and numerical solutions (eqs. 4 and 6) obtained 
using trigonometric relations is e(k) = 2 sin{k[a — 
a(k)]t / 2). Expanding the error in small phase angles 
yields the simple expression for the phase error: 


where T 0 is the initial time, and C and L are problem 
dependent constants related to solution smoothness, 
etc. In equation (2) the time-step satisfies At < 1, 
and a given error tolerance can be achieved by in- 
creasing the order p while increasing the time-step. 
The resulting algorithm is more efficient if any addi- 
tional work accrued at each large time-step, is more 
than compensated by a reduced number of steps. 
Work, however, increases with the order p. Near op- 
timal values of p in the range 3 < p < 5 exist 
for countless stiff and non-stiff model problems, (see 
§IV. 10 in Hairer and Wanner 56 ). 

While it is generally recognized that high-order 
temporal schemes result in greater efficiency for 
time-dependent, problems, it is less well appreciated 
that higher order spatial schemes additionally con- 
tribute to time-dependent, efficiency! The virtues of 
high-order spatial schemes were first recognized and 
quantified by Kreiss and Oliger. 83 To illustrate this 
advantage we provide an overview of the original ar- 
gument presented in Kreiss and Oliger. 83 

Consider the wave equation and initial data 


e(k) = | k[a - a(k)]t | (9) 

Taylor series arguments produce the leading order 
term for the difference in wave speeds 

a — a(k) ~ ad p {kAx) p (10) 

with 0 P a scheme and order dependent constant. 
Substituting equation (10) into (9), defining the 
points per wavelength as P = N/k, and substi- 
tuting t = Tf yields 

e(k) ~ [,3 p kaT f (^) P ] (11) 

The semi-discret.e solution error accumulates lin- 
early with time Tf, and is a strong function of the 
spatial truncation error. Rearranging equation (11) 
in terms of a maximum acceptable target error exile) 
yields the expression: 


Pp > 2 7T 



( 12 ) 


Ut + aU x = 0, U (x, 0) = e ikx (3) 

on the space and time intervals 0 < * < 2 tt and 
0 < < < T) , with the exact solution: 

U{x,i) = e ik ( 4 ) 


The grid points per wavelength necessary to achieve 
the specified target error ex (k) increases for prob- 
lem size (2 tt). The other dependencies rapidly de- 
crease as the order of the spatial approximation is in- 
creased. motivating high-order spatial formulations 
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in time-dependent problems. The cost of the com- 
putations increase with increasing order of accuracy, 
and a global minimum is reached for a finite value 
of p. Kreiss and Oliger 83 suggest 4 < p < 6 for 
problems of practical interest. Note that the optimal 
order for spatial and temporal operators is similar. 

In summary, error accumulates linearly in time. 
The global error at Tf is the sum of local errors 
that accumulate from each time-step in the integra- 
tion. The local error at each time-step is the sum 
of three components: the temporal truncation error, 
the spatial truncation error, and the algebraic er- 
ror. A simulation requiring many time-steps to reach 
Tf requires extremely small local errors. High-order 
methods are the most efficient means of achieving 
these small local error tolerances. 

An example will help to clarify this point. Con- 
sider a steady-state problem requiring lift to an engi- 
neering accuracy of three significant digits. A second 
order method could easily achieve this accuracy re- 
quirement. Now consider a similar time-dependent 
problem requiring lift (at the specified Tf) to an en- 
gineering accuracy of three significant digits. Fur- 
ther assume that 100 time-steps are required to inte- 
grate from To to Tf . The local error (temporal, spa- 
tial, algebraic) at each time-step must be less than 
10 -r ’ to achieve the desired error tolerance of three 
significant digits. The constraint ou spatial error 
(10 -5 ) in the time-dependent problem is much more 
severe than that required for the steady-state prob 
lem (10 -3 ). This example demonstrates the com- 
pelling need for high-order spatial operators espe- 
cially for time-dependent simulations. 

3.2 Temporal A l gorithms 

3.2.1 Overview 

The application of method of lines (MOL) to time- 
dependent partial differential equations (PDEs) re- 
sults in an initial value problem (IVP) for a system 
of ODEs. Dozens of excellent texts with detailed de- 
scriptions of multi step, multi stage, and linear multi 
step methods have been written on the numerical 
integration of ODEs. 26 - 40 - 55 - 56, 86, 128 After more 
than 100 years of theoretical development, the math- 
ematical framework for solving ODEs is relatively 
mature. In a general context, it is doubtful that 
dramatic (factors of 10) efficiency improvements can 
come from new methods. 

The potential for dramatic efficiency 7 improve- 
ments is greater in the field of time-dependent. CFD, 
where current methodologies are surprisingly primi- 
tive. This schism between tidy mathematical theory, 
and rough CFD practices is not without good reason. 
Fluids practitioners are preoccupied with more ur- 


gent issues such as algebraic solvers, dimensionality 
issues, discontinuities, nonlinear instability, turbu- 
lence models, grid generation, etc. Nevertheless, the 
current objective is to identify mature technologies 
in the ODE literature that could have an immediate 
impact in CFD. 

The hallmark of current ODE software is the abil- 
ity to perform automated integration for stiff ODEs. 
The first widely available multi step integration li- 
brary was that developed by Gear, 51 later modi- 
fied and improved by Hindmarsh, 60 resulting in the 
LSODE family of codes. Other variants have prolif- 
erated over the past two decades to account for the 
deficiencies of the original approaches (see VODE 
61) - 

Automated integration begins with the user spec- 
ifying 1) the system of ODEs, 2) the system Ja- 
cobian, and 3) the desired solution error tolerance. 
The software then automatically integrates the equa- 
tions, using the most efficient numerical method cho- 
sen from a variety of candidate methods (second- 
order backward differentiation formulae (BDF2) is 
frequently used). A reliable solution error estima- 
tor allows variable time-stepping. The time-step is 
adjusted to match the desired error tolerance. The 
resultant nonlinear system of algebraic equations is 
solved at each time-step using a Newton or modified 
Newton method. Direct matrix inversions are used 
within the Newton methods whenever possible. The 
algebraic error is reduced to a predetermined level, 
a constant multiple below the specified error toler- 
ance. The Jacobian used in the nonlinear iteration 
is periodically reevaluated and stored based on the 
convergence rate of the iteration. 

In contrast, the second-order accurate multi step 
BDF2 method is extensively used in the CFD com- 
munity. System dimensionality prohibits the use of 
direct inverse methods useful for Newton or modi- 
fied Newton methods. Iterative techniques such as 
Newt.on-Krvlov methods are usually not as efficient 
as other more highly 7 tuned methods (multigrid or 
combinations of methods) . Error estimation or vari- 
able time-stepping mode is not perceived as neces- 
sary 7 (correct or otherwise). 

Three technologies presently used in the ODE 
community could have an immediate impact on 
CFD: 

• high-order integrators : p > 3 

• error estimation/variable time-stepping 

• iteration termination strategies 

To support this assertion, a brief summary 7 of each 
area is presented. 
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3.2.2 High-Order Integrators 

All general purpose solvers must integrate equa- 
tions of considerable stiffness. We begin with a 
broad overview of stiffness, and identify the mathe- 
matical properties that enable a temporal integrator 
to efficiently integrate stiff equations. 

Consider the integration of the system of ordinary 
differential equations represented by the equation 


f = S(U( ^ 


In the present case, the vector S results from the 
semi-discretization (spatial and source terms) of the 
equations of fluid mechanics plus a suitable turbu- 
lence model. The integrator must integrate any S 
with which it is provided. Numerical difficulties of- 
ten arise when the Jacobian of S, ./ = dS/d\J, has 
large eigenvalues. A useful definition for stiffness 
states that a problem is stiff when the largest eigen- 
value of the Jacobian J (scaled by the time-step) 
\\z = A(Af)|| contained in the complex left-half- 
plane (LHP) becomes much greater than unity. The 
resulting stiffness is then governed by both the Jaco- 
bian and the chosen time-step. Ideally, the time-step 
is selected solely based on error considerations and 
a good method simply executes this step-size in a 
stable and robust fashion. Time integration meth- 
ods that do not amplify any LHP scaled eigenvalues 
are called A-stable. While A-stability is generally 
necessary, it is often not sufficient. We further de- 
mand that all eigenvalues ||« — > — oo|| be completely 
damped. The combination of these two properties, 
A-stability and damping of — oo eigenvalues, is de- 
fined as L-stability. General purpose solvers invari- 
ably rely on L-stable methods (and the partially sta- 
ble L{a) methods with suitable error controllers) to 
suppress temporal numerical instability and facili- 
tate convergence of the nonlinear equation solver. 

Popular implicit ODE integration methods are 
generally either distinctly multi step or multistage 
methods. Each has different strengths and weak- 
nesses. Implicit multi step BDF methods compute 
each U-vector update to design order of accuracy 
using one nonlinear equation solve per step. Unfor- 
tunately, they are not A-stable above second-order. 
Additionally, they are not self-starting and have di- 
minished stability properties when used in a vari- 
able step-size context. (Stability proofs are formu- 
lated assuming constant time-steps. Variable time- 
step cases may not be stable.) Practical experience 
indicates that large-scale engineering computations 
are seldom stable if run with BDF4. 102 The BDF3 
scheme, with its smaller regions of instability, is 
often stable but diverges for certain problems and 


some spatial operators. Thus, a conservative prac- 
titioner uses the BDF2 scheme exclusively for large 
scale computations due to its L-stability rather than 
L(a)-stability. 

Practical Runge-Kutta (RK) methods such as ex- 
plicit, singly diagonal implicit, Runge-Kutta (ES- 
DIRK) methods can be made arbitrarily high-order 
while retaining L-stability but possess intermedi- 
ate U- vectors wi III a reduced order of accuracy and 
lesser stability. This reduced stage order may give 
rise to order reduction phenomena in the presence 
of substantial stiffness. ESDIRK schemes with s 
stages require (s-1) nonlinear equation solves per 
step. Achieving progressively higher stage-order 
methods is possible with fully implicit methods such 
as the Radau IIA family. The cost of fully implicit 
methods greatly exceeds that of ESDIRK methods 
in the current context. Much less experience exists 
with implicit RK methods than BDF methods in the 
computation of large-scale engineering flows. 

The general formula for a k- step, order-A, BDF 
scheme can be written as 
fc-i 

U (n+fc) = _^ aiU (»+0 + (A<)/?fcS^ +fc (>13) 

where n is the time-step index. At each time-step 
the BDF involve the storage of k+ 1 levels of the so- 
lution vector U, and the implicit solution of one set 
of nonlinear equations. Stability diagrams for these 
methods may be found in Hairer and Wanner. 56 At 
order k > 2 an unstable zone for scaled eigenvalues 
in the complex LHP exists. At orders {1,2, 3, 4, 5, 6} 
the methods are L(a)-stable where a is given by 
{90°, 90°, 86.03°, 73.35°, 51.84°, 17.84°}. (The dif- 
fusion terms yield negative real eigenvalues. The 
convective terms discretized with nondissipative op- 
erators yield eigenvalues clustered on the imaginary 
axis a = 90°. Numerical dissipation and boundary 
conditions displace both sets of eigenvalues into the 
LHP. Spatial operators with high levels of dissipa- 
tion are more likely to be stable with BDF3.) 

ESDIRK methods 78 ' 85 are implemented as 

k 

U fc = U" + (At) ^ a k jS (U J ) , k = 1 , s 
l=i 

U ra+ i = U" + (At)'^Tb j S(U j ) (14) 

i=i 

U a+1 = U" + (W) 

i=i 

where s is the number of stages, a^j are the stage 
weights, bi and bj are the main and embedded 
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scheme weights. The vectors U and U are the 
p tft -order and ( p — l) tA -order solutions at time level 
n + 1. The vector U is used solely for estimat- 
ing error and is calculated at little extra cost. ES- 
DIRK schemes differ from traditional SDIRK meth- 
ods (see §IV.6 in Hairer and Wanner 56 ) by the 
choice an = 0, which permits stage-order two meth- 
ods. The stiffly accurate assumption ( a s j = b;) 
makes the new solution U n + 1 independent of any 
explicit process within the integration step. 

Many methods that combine multi step and multi 
stage schemes, generally referred to as linear multi 
step (LMS) methods, have been proposed in attempt 
to overcome the problems of each. Surprisingly, the 
results obtained with LMS methods (with a few ex- 
ceptions) have been disappointing compared with ei- 
ther multi step or multi stage schemes. Notable LMS 
methods include the works of Cash. 34 ' 35 To increase 
the stability of the BDF methods, Cash proposed the 
extended backward differentiation formulae (EBDF) 
and the modified EBDF (MEBDF) schemes. The 
MEBDF schemes involve three stages to advance the 
solution one time-step. The first two stages are built 
from existing ( p — l) t/i -order BDF formulas, while 
the last stage combines the two previous BDF results 
into ap t/l -order solution. Note that the second BDF 
stage predicts a “super-future” point one time-step 
beyond the target time level, and substantially con- 
tributes to the A-stability of the method. At orders 
{1, 2, 3, 4. 5, 6} the methods are L(a)-stable where a 
is given by {90°, 90°, 90°, 90°, 88.36°, 83.07°}. The 
machinery involved with implementing the MEBDF 
algorithm is nearly identical to that involved in the 
BDF formulations. MEBDF schemes have the ad- 
ditional advantage that very accurate solution data 
are available on the first and third stages, based on 
previous information. This information can be used 
to provide the starting guess for the nonlinear it- 
eration, and to establish time-step error estimates. 
The second stage typically uses the trivial guess as 
the starting point for the nonlinear iteration and no 
error estimate is made. 

Butcher 28 proposed a class of LMS methods for 
stiff differential equations. The new methods com- 
bine the properties of A- and L-stabilit.y, and are rea- 
sonably simple to implement. They have a stability 
region that, is identical to that of a R.K method, but 
have high stage order. Uniformly high stage order 
eliminates the possibility of order reduction. The 
new methods were identified by focusing specifically 
on a diagonally implicit subclass of schemes referred 
to as DIMSIm". 27 


3.2.3 Error Estimation 

Temporal error management in the CFD com- 
munity is presently accomplished by systematically 
halving the time-step until the solution is indepen- 
dent of further reduction. This strategy, while ac- 
complishing the desired goal, can be streamlined by 
using an error estimator at each time-step and ad- 
justing each time-step to attain the desired error. 

Error estimation is accomplished by comparing 
two solutions of different orders (U" + 1 and U n+1 ) 
at the same time-step. For reasons of efficiency, the 
auxiliary solution U n+1 should be available at little 
additional cost. For example, in ESDIRK schemes 
(see eqn. 15), as well as MEBDF 35 schemes, both 
U n+1 and U n+1 are constructed from available data. 
The difference ||U’ I + 1 — U ra + 1 1| is proportional to the 
truncation error of the lower order formula U ,!+1 . 
The estimate predicts the magnitude of the error in 
the solution, and gives insight into its overall qual- 
ity. Frequently, linear and nonlinear instability can 
be predicted by the estimator well before the simu- 
lation diverges. 

Figure (1) show's the error estimate (MEBDF4) 
for various At. The test problem is for periodic 
shedding from the turbulent circular cylinder. The 
estimates are accurate to the correct order based 
on grid-converged data. The error estimate predicts 
that certain portions of the shedding cycle are more 
difficult to resolve in time. Variable time-stepping 
could easily increase the efficiency of the calculation 
by adjusting the time-step so that the same amount 
of error is produced at each time-step. 

Variable time-stepping can introduce instability 
into some temporal integrators. The stability func- 
tion of multi step schemes (BDF, MEBDF, LMS) 
is derived assuming constant time-steps. Large de- 
partures from constant step-size can lead to solution 
instability (although a good error estimator should 
forewarn this possibility). Conversely, the stability 
of multi stage schemes (ESDIRK) is independent of 
variable time-steps, because they are self-starting. 

The time-steps in a variable time-step formula- 
tion, are chosen by a controller. A simple explicit 
controller is (see §IV.2 in Hairer and Wanner 5S ) 


(A t ) n+1 = (A t) n (- 


Tol 

U"+! — U n+1 | 


(15) 


Similar yet more elaborate controllers exist for im- 
plicit formulations. The stability characteristics of 
a controller can be tuned/optimized in conjunc- 
tion with the integration technique it is control- 
ling. Together, they should meet the design ob- 
jective and not introduce instability into the inte- 
gration. Note that in figure (1) the predicted error 
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at finer tolerances has a high frequency component, 
that the controller must suppress. (See Kennedy 
and Carpenter' 8 for details on the feedback error 
controllers used with the ESDIRK scheme.) 

3.2.4 Termination Strategy 

An accurate error estimate can also be used to 
automate the termination strategy of the nonlinear 
iteration. Two competing components of temporal 
error are the truncation and algebraic errors. Trun- 
cation error is related to At and the order of ac- 
curacy p, while algebraic error is the residual error 
generated each time-step by approximately solving 
the algebraic system. The local temporal error is 
the sum of the two components. To see full design 
order from the temporal scheme, the algebraic er- 
ror must be driven below the truncation error at 
each time-step. This requires an accurate measure 
of truncation error, and must be provided by the 
error estimator. 

The iteration termination strategy is complicated. 
Our experience indicates that design-order temporal 
convergence is achieved by maintaining a tolerance 
ratio of 10 -J < T < 10“ 4 . Here T is defined as the 
ratio of nonlinear algebraic error to temporal inte- 
gration error at each time-step (or stage). Algebraic 
error for the nonlinear iteration is based on the L 
norm of the density residual. Choosing the time-step 
based on accuracy considerations alone may not, be 
the most efficient strategy for a temporal calcula- 
tion. Decreasing the time-step can possibly greatly 
increase the convergence rate of the nonlinear alge- 
braic system, thus increase efficiency. Gustafsson 
and Soderlind 54 devised optimal criteria for adjust- 
ing At. They assumed that either fixed point, itera- 
tions, or modified Newton iterations is used for solv- 
ing the algebraic system. The time-step is adjusted 
so that the iteration convergence rate approximately 
equals the optimal value. Because typical CFD alge- 
braic solvers fall somewhere between fixed point, and 
modified Newton iterations, additional work to re- 
fine these estimates is needed in the context of CFD 
time-dependent, solvers. 

3.2.5 Bottlenecks 

Algebraic solvers that exhibit, poor convergence 
behavior are an impediment for high-order schemes. 
Huge time-steps are needed to utilize the favor- 
able aspects of high-order formulations. Algebraic 
solvers are needed that, exhibits time-step indepen- 
dent convergence characteristics. If the convergence 
rate varies considerably with the time-step, then it, 
may be more efficient to use a low-order scheme 
with small time-steps. Thus, high-order temporal 
schemes need fast and robust, algebraic solvers. Tur- 



Figure 1. Time-dependent variation of predicted 
Lo density error as calculated with the MEBDF f 
scheme. The test case is turbulent flow around a 
circular cylinder, at Re = 10 4 , and Ma = 0.25. 

bulent cases that have little or no convergence (one 
order of magnitude) present, a second obstacle. A 
small number of cases are extremely difficult to con- 
verge, yielding dubious solutions at best. Neverthe- 
less, solutions are still sought,. It, may be difficult 
to keep high-order formulations stable under these 
circumstances. 

A perceptual impediment is the implementation 
of error estimation technology. A change in atti- 
tude about the nature of temporal error and the 
importance of its control is necessary. In spite of 
the perceived adequacy of existing temporal error 
practices, the CFD community should immediately 
adopt the practice of reporting a time-step error es- 
timate as a necessary requirement of a high fidelity 
time-dependent simulation. Ideally, the estimate 
should include the component of primary interest 
in the simulation. For example, if lift and drag are 
the object of the study, then the estimate should in- 
clude stepwise error estimates of these quantities, as 
well as information on which formulas were used to 
obtain the estimate. 

Another common yet dangerous practice in the 
CFD community is to use a fixed number of itera- 
tions for each time-step. This approach eliminates 
the need for an iteration termination strategy and 
in most circumstances is satisfactory. Global tem- 
poral error [see equation (2)] strongly depends on 
time-steps with large local error. The error from just 


6 



one nonconvergent time-step potentially could domi- 
nate the error from all other time-steps combined! If 
an intrinsic feature of the flow significantly changes 
the convergence rate of the algebraic solver, then a. 
fixed number of iterations is not a good strategy. 
The periodic blowing from zero mass fluid actuators 
is a prime example. Different phases of the cycle 
converge at different rates because convergence rate 
is sensitive to boundary conditions. A termination 
strategy that ensures a uniformly bounded algebraic 
error at each time-step is needed. 

3.2.6 Langley effort 

An ongoing effort focuses on the efficacy and effi- 
ciencies of several time integration schemes for the 
unsteady compressible Navier-Stokes equations. Ex- 
isting and newly developed multi step and multi 
stage schemes are being studied, with particular at- 
tention to high-order (p > 3) schemes. Past work 
includes comparisons of the high-order (ESDIRK4) 
7 8 Runge-Kutta scheme with first- and second-order 
BDF on laminar problems. Bijl, et al. 20 showed that 
the efficiency of the ESDIRK4 scheme exceeds that 
of the BDF2 by a factor of 2.5 at engineering er- 
ror tolerance levels (10 _1 -10 -2 ). Efficiency gains are 
more dramatic at smaller tolerances. No problems of 
nonlinear instability were noted with the high-order 
ESDIRK4 scheme on the problems tested. 

Carpenter et al. 32 has shown that stage order two 
Runge-Kutta schemes are susceptible to order re- 
duction for stiff systems, although none is experi- 
enced for laminar problems with stiffness levels of 
0( 10 3 ). However, turbulence models exhibit con- 
siderable stiffness at Reynolds numbers in the range 
of 10 5 — 10'. Significant order reduction is expe- 
rienced with ESDIRK4 for cases experiencing stiff- 
ness from strong turbulence fields. Ongoing stud- 
ies include investigating the efficiency of ESDIRK 
schemes on other one- and two-equation turbulence 
models. 

Figure 2 shows the convergence behavior of the 
ESDIRK4 scheme, the BDF2 and BDF3 schemes, 
and the MEBDF4 scheme. The test problem is the 
circular cylinder at Reynolds number 10 4 , with a 
Mach number of 0.25. The calculations are run 
with the unstructured Fun2D code. 2 Design order 
slopes are obtained for each scheme: 2, 3, 4, and 
3 for BDF2, BDF3, MEBDF4, and ESDIRK4, re- 
spectively (note that w'e have accounted for the the- 
oretical order of ESDIRK4 in accordance with order 
reduction). The MEBDF4 scheme was added to the 
comparison because it is a stage order three method 
and is not as susceptible to order reduction as the 
ESDIRK4 scheme. Research continues on establish- 
ing reliable error estimators and iteration termina- 



Figure 2. A comparison of temporal density error 
obtained with BDF2, BDF3, MEBDFf. and ES- 
DIRK \ schemes on a circular cylinder, with Re 
= 10 4 , and Ma = 0.25. The turbulence model is 
Spala rt- Alima ras. 

tion strategies. A comparison of efficiency will be 
made between all integration schemes once each is 
automated . 

A final observation is relevant to help focus future 
work. The BDF2 scheme yields engineering accuracy 
if each temporal mode in a time periodic flow is re- 
solved with approximately 50 — 100 time-steps. The 
fourth-order ESDIRK formulation attains a similar 
accuracy (in spite of order reduction) using 5 — 10 
time steps per period, with five stages per step, 
yielding approximately 25 — 50 time-samples/period. 
State-of-the-art LMS methods could lower this esti- 
mate to 15 — 30 time-samples/period, an improve- 
ment of approximately 0( 10 1 / 2 ) over existing tem- 
poral efficiency. The theoretical lower bound for 
temporal schemes, based on infinite order Chebyshev 
operators, is 7r samples/period. High-order schemes 
are presently asymptotically close to this theoreti- 
cal lower bound. Another factor of three reduction 
in samples/period is perhaps all that remains and is 
becoming increasingly more difficult to attain. Al- 
gorithmic work focusing on other aspects of solver 
technology will have a greater chance of producing 
meaningful improvements during the next decade. 
The next section describes high-order spatial algo- 
rithms and their potential to increase the temporal 
efficiency, and section (4) describes the current and 
future status of algebraic solvers. 
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3.3 Spatial Algorithms 

3.3.1 Overview 

The spatial algorithms used currently in general 
purpose aerodynamics solvers have not changed ap- 
preciably during the past decade. Most current pro- 
duction codes (structured or unstructured) rely on 
some form of second-order upwind formulation with 
flux limiting to provide necessary robustness in the 
vicinity of unresolved features in the flow. Excellent 
texts describing these methodologies can be found 
elsewhere. For finite-difference methodologies see 
Hirsch 62 ’ 63 , and LeVeque. 88 For basic finite-element 
methodologies see Hughes 70 , Zienkiewicz and Tay- 
lor 149 ' 150 and Baker and Pepper. 10 

In section 3.1, we established that spatial algo- 
rithms play a important role in determining tem- 
poral efficiency. High-order methodologies will sig- 
nificantly contribute to the ultimate goal of effi- 
cient, general-purpose, time-dependent aerodynamic 
solvers. A broad overview of the spatial discretiza- 
tion landscape is now presented. Enabling technolo- 
gies that allow extension of high-order methods into 
the general purpose aerodynamic solver arena are 
identified. A wealth of scientific literature supports 
the assertion that high-order general purpose algo- 
rithms will likely mature within the unstructured 
finite element framework within the next decade. 

Spatial operators are categorized by the kind of 
grids on which they are formulated: structured 

or unstructured. A structured grid has large re- 
gions of the interior vertices that are topologically 
alike, which results in well-established connectiv- 
ity patterns. Accommodation of complex geome- 
tries requires an arbitrary subdivision of the struc- 
tured grid into what is referred to as a hybrid or 
multi block formulation. Three examples of struc- 
tured codes currently used at Langley as general- 
purpose aerodynamics solvers include the block- 
structured TLNS3D. 144 and CFL3D, 142 and the 
overset-structured OVERFLOW 72 . An unstruc- 
tured mesh is one in which vertices may have ar- 
bitrarily varying local neighbors. Three examples 
of unstructured codes used at Langley (ICASE) are 
USM3D 49 , FUN3D 2 , and NSU3D. 94 The distinction 
between structured and unstructured meshes usu- 
ally (although not necessarily) extends to the shape 
of the elements: 2-D structured meshes typically use 
quadrilaterals, while unstructured meshes use trian- 
gles, with similar analogous element shapes in 3-D 
(hexahedra vs. tetrahedra). 

Structured solvers offer simplicity, easy data ac- 
cess, and thus efficiency. The data structure and 
algorithmic simplicity of structured solvers leads to 
more efficiency and lower memory requirements for 


a given accuracy tolerance. A discrete derivative 
requires simple increments/decrements in array in- 
dices, in stark contrast to an unstructured formu- 
lation. The structured advantage in CPU time and 
memory can be as much as a factor of three on prob- 
lems not requiring significant grid adaptation. How- 
ever, on a complicated geometric domain a struc- 
tured mesh may require many more elements than 
an unstructured mesh, because elements in a struc- 
tured mesh cannot vary in size as rapidly. Struc- 
tured grid generation approaches are far from being 
fully automated, and require user guidance in the 
decomposition step. A complicated 3-D structured 
mesh can take a month to generate. The current and 
future role of structured formulations is for repeti- 
tive computations late in the design cycle where grid 
templates might exist and grid generation and adap- 
tation are not important components in the solution 
process. 

Unstructured meshes offer flexibility in fitting 
complicated domains, rapid variation from small 
to large elements, and relative ease in refinement 
and de- refinement. Unlike structured mesh gener- 
ation, unstructured mesh generation has been au- 
tomated in mainstream computational geometry for 
some years. The major approaches for generating, 
refining, and improving unstructured meshes rely 
on unconstrained and constrained Delaunay trian- 
gulation, quad trees algorithms, or combinations of 
the above. 19 A highly effective combination of tech- 
niques for high Reynolds number flows is an advanc- 
ing layers method (ALM) 114 in the near- wall re- 
gion, and an advancing front method (AFM) 91 in 
the far- field. Highly stretched viscous grids can be 
generated in a reasonably automated fashion with 
this approach. Automation begins to break down as 
aspect ratio increases on complex geometries. 

Element shape has a. profound impact on the accu- 
racy and efficiency (direct and indirect) of a formu- 
lation. Meshes with unintended large aspect ratio 
cells lead to both poorly conditioned matrices and 
poor solution accuracy. Poor solution accuracy re- 
quires more grid points for a given accuracy. The 
additional cost of fixing a bad mesh can usually 
be mitigated by the faster convergence of the iter- 
ative solver. 1, Babuska and Aziz 7 showed that con- 
vergence on triangular elements is achieved only for 
angles bounded away from 180°. This rather weak 
condition becomes an issue for strongly anisotropic 
meshes used in high Reynolds number turbulent 
Navier-Stokes simulations. Near-wall aspect ratios 
on these grids can be in the range 10 4 — 10 5 . For- 
mulations typically try to limit the maximum an- 
gle in a grid (for example 179°) even though cur- 



rent evidence is divided on the necessity of this con- 
dition. Quadrilateral and hexahedral meshes have 
an advantage in accuracy over triangular and tetra- 
hedral meshes for these problems. The faces of 
hexahedral elements in the boundary layer are ei- 
ther almost parallel or almost orthogonal to the sur- 
face. Shock fronts and shear layers, which are also 
strongly anisotropic, require high aspect ratio cells 
for which the direction and location cannot be pre- 
dicted in advance. Generation of these meshes can 
be difficult. 

General purpose aerodynamic solvers have pro- 
gressively shifted from hybrid/structured methods 
to unstructured formulations over the past decade. 
The principal motivations driving this change are 
grid generation on complex configurations and grid 
adaptation. These compelling reasons are likely to 
become more important during the next decade, par- 
ticularly as the grid adaptation field matures for 
time-dependent simulations. 

A large amount of inertia persists in the struc- 
tured grid world, which is not entirely counterpro- 
ductive. A time-dependent niche exists for computa- 
tionally efficient formulations over the next decade. 
Unlike 3-D steady-state computations, realistic 3- 
D time-dependent computations are presently con- 
strained by processor speed rather than memory re- 
quirements. The increased efficiency of structured 
methods is a notable advantage when run times can 
be decreased by a factor of two to three. The addi- 
tional hybrid/structured grid generation time can be 
amortized if a calculation is likely to run for months. 

3.3.2 High-Order Spatial Operators 

High-order spatial operators need fewer points 
than second-order operators, to resolve the same in- 
formation. The exact reduction strongly depends 
on the desired accuracy. Steady-state problems re- 
quiring an accuracy of three significant digits can be 
achieved with fourth-order schemes in half as many 
points in each spatial dimension. The total reduc- 
tion in the number of points is approximately 0( 10 1 ) 
in 3-D. Time dependent simulations that require so- 
lution accuracy to four or even five significant digits 
at each time-step, will favor high-order formulations 
to a larger degree. High-order spatial methods can 
increase the efficiency the time-dependent simula- 
tions by 0( 10 1 - 10 2 ). 

The constraints necessary to expedite grid gener- 
ation and grid adaption will guide the next genera- 
tion of high-order solvers. High-order methods must 
move beyond proof of concept and into the realm of 
being tools used to increase the efficiency of aerody- 
namic solvers. 


The implementation of high-order methods is 
strongly dependent on whether the grid is structured 
or unstructured. High-order finite-elements (FE) are 
natural candidates for structured or unstructured 
meshes, while high-order finite-difference (FD) tech- 
niques are usually implemented on block-structured 
or overset grids. Finite- volume (FV) techniques ex- 
ist in both forms; a close similarity between linear 
element FE methods and FV methods exists. All 
three approaches solve different, forms of the govern- 
ing integral equation. FV directly solves the integral 
equation by approximating the numerical fluxes. FD 
solves the divergence form of the integral equation 
by approximating the derivatives. FE take the di- 
vergence of the integral equations, multiply by an 
arbitrary test function, and integrate by parts. The 
solution itself is the resulting approximation. 

Not all current general purpose spatial discretiza- 
tion algorithms are natural candidates for high- 
order extensions. For example, based on 2-D re- 
sults. Casper and Atkins 36 noted that a 3-D hy- 
brid/structured essentially nonoscillatory (ENO)- 
FV formulation would be extremely expensive to 
implement relative to comparable FD techniques. 
Barth and Frederickson 12 , and Barth 13 extended 
their unstructured FV solver to account for A;-exact 
reconstruction. They note that comparing quadratic 
with linear reconstruction on a triangle requires 
roughly quadruple the number of solution unknowns. 
All high-order formulations require more work than 
second-order formulations, but some are more effi- 
cient than others. 

Many different approaches to high-order FE have 
been adopted in developing numerical schemes to 
solve the compressible Euler equations. Two major 
classes have emerged as candidate schemes: 1) sta- 
bilized methods (continuous across interfaces), and 
2) discontinuous methods (discontinuous across in- 
terfaces) . 

Standard Galerkin FE dis- 

cretizations of convection-dominated Navier-Stokes 
equations produce wildly oscillating solutions un- 
less dissipation terms are added to the formulation. 
Since the early 1980s stabilized FE methods have 
become increasingly popular in CFD. Early devel- 
opment motivated by the success of upwind FD/FV 
schemes included the streamline diffusion finite el- 
ement Method (SDFEM), 68 ' 75 which later evolved 
into the streamline upwind Petrov-Galerkin (SUPG) 
scheme of Brooks and Hughes. 25 A stabilizing term 
is added into the weak statement motivated by in- 
viscid terms and results in a perturbed standard 
Galerkin test function. The stabilization creates 
an upwind effect by weighting more heavily the up- 
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stream nodes within each element. Hughes and Tez- 
duyar 69 generalized the SUPG method to first-order 
hyperbolic systems and included work on the Euler 
equation. The original SUPG formulation suffered 
from oscillation in steep gradient regions (shocks) 
which led to the introduction of entropy variables 
and ultimately to stabilization terms including the 
effects of both the inviscid and viscous terms in their 
Galerkin least-squares (GLS) formulation. 71 A vari- 
ety of methodologies have been proposed to provide 
additional stability to the convection terms, mono- 
tone discrete solutions and ease of implementation. 

Another approach that has been gaining in pop- 
ularity in recent years is the discontinuous Galerkin 
(DG) method. The DG method originally intro- 
duced by Reed and Hill 117 exhibits several distinct 
advantages when applied to complex unstructured 
grids. Local polynomials are used to represent the 
data to arbitrary order, with the data on element 
interfaces treated as discontinuities. The approach 
is advantageous because the solution accuracy is rel- 
atively insensitive to mesh smoothness and can be 
extended to arbitrarily shaped elements. In 1986, 
Johnson and Pitkarata 76 proved that the conver- 
gence rate of the method is ( Ax) k+l ^ for general tri- 
angulations. The method generates a local mass ma- 
trix that can easily be inverted, making the method 
efficient for explicit time integration. An entropy 
inequality for any scalar nonlinear equation 73 ex 
ists, implying discrete nonlinear instability for dis- 
continuous solutions. (This assumes wellposedness 
and boundedness of the continuous nonlinear prob 
lem.) Several researchers have demonstrated super- 
convergence with DG. 37 Lowrie et, al. 92 obtained 
convergence rates of 2 p + 1, and Hu and Atkins 67 
showed that the dispersion of the DG method is gov- 
erned by 2p T 1 for polynomials of order p. The DG 
formulation produces a matrix with many dense but 
small sub-matrices weakly coupled to their neigh- 
bors. Algorithmically, this matrix structure excels 
in a parallel environment and has high cache effi- 
ciency. Atkins and Shu 6 developed a quadrature- 
free approach that allows the precomputation and 
storage of much of the algorithm, thereby increasing 
the efficiency. 

In the early 1980s, the p- and the hp-FEM meth- 
ods were introduced by Babuska, and Szabo 8 . They 
showed that for elliptic problems exponential conver- 
gence could be achieved with the hp-FEM method. 
The degree of the approximating polynomial can 
vary by elements so both grid refinement and order 
refinement are used simultaneously to attack solu- 
tion error. The methods show design order for p 
fixed in the limit h — > 0, and convergence for h 


fixed p — > oo. The behavior for high Reynolds num- 
ber Navier-Stokes equations is less clear; neverthe- 
less the hp-FEM methods have great potential in the 
context of complex geometries and grid adaptation. 

High-order FD excel in their simplicity and effi- 
ciency, and in the richness of linear and nonlinear 
algorithmic permutations that can be formulated. 
The Achilles heels of FD are “boundaries” and “sta- 
bility” . 

Achieving numerical stability near boundaries 
with high-order FD stencils is difficult. This insta- 
bility is closely related to the classical Runge oscil- 
lations exhibited by high-order polynomials on uni- 
form grids near the boundaries. The solution in both 
cases is to lower the polynomial order, compress the 
grid, or increase the stencil width. Gustafsson 53 
showed that to maintain spatial design order accu- 
racy, the boundary stencil order must not deviate 
by more than one from the interior order of accu- 
racy: fourth-order interior stencils require boundary 
closures of at least third-order accuracy! Note that 
in fluid dynamics applications, near- wall regions are 
precisely the regions were high-order accuracy is de- 
sirable. In general, low order treatments for bound- 
ary closures are not an acceptable alternative unless 
special near-wall grid refinement is used to compen- 
sate for reduced accuracy. 

Strand 137 following the work of Kreiss and 
Scherer. 84 partially resolved the boundary closure 
dilemma by presenting constructive procedures for 
developing stable and accurate boundary schemes. 
Stability is ensured in an L 2 norm using a dis- 
crete summation-by-parts (SBP) procedure. Car- 
penter et al 30 and Olsson 109, 110 showed how to 
impose the physical boundary conditions to pre- 
serve the SBP energy estimate. To date, bound- 
ary closures have been formulated for central and 
upwind FD schemes and Hermitian compact FD. 
FD schemes require smooth structured meshes that 
are often difficult to generate on complex geome- 
tries. Multi block/overset grids relax the gridding 
constraints, allowing piecewise smooth grids around 
complex geometries, but create a new set of compli- 
cations. Conservation is a major concern on multi 
block grids and is extremely difficult to achieve on 
overset grids. Shocks and other discontinuities near 
interfaces must be treated carefully to ensure correct 
shock speeds and locations. In addition, interfaces, 
like boundaries, can cause linear and nonlinear in- 
stability and lead to decreased levels of solver ro- 
bustness. 

Multiple attempts have been made to overcome 
the difficulties of complex geometries for high-order 
FD schemes. By far, the most common solution to 
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ensure conservation and stability has been to reduce 
boundary or interface accuracy. As a general rule, 
this approach works well if little structure exists near 
the boundary or interface. Solution accuracy in com- 
plicated flow scenarios, however, is difficult to pre- 
dict. Carpenter, Nordstrom, and Gottlieb 31 > 106, 107 
have developed instable interface conditions based 
on SBP energy estimates. The interface points are 
treated discontinuouslv through a penalty term and 
provide conservative, high-order solutions on multi 
block grids. The only grid requirement is Co inter- 
face continuity between blocks, a mild restriction. 

3.3.3 Bottlenecks 

The major obstacle facing all high-order spatial 
discretization methods FE or FD, structured or un- 
structured, is nonlinear instability in the presence 
of unresolved features. Shock and sliplines are no- 
table examples. The classical approaches to deal 
with discontinuities in FD and FV are the addition 
of local artificial viscosity and/or filtering, and total 
variation diminishing (TVD) or limiting approaches. 
Equivalent approaches exist in FE, though they are 
termed “stabilization.” The amount of added dissi- 
pation depends on the simulation objectives. Mono- 
tone solutions can be obtained with any formulation 
at the expense of reduced accuracy. In principle, the 
minimum amount of dissipation necessary for non- 
linear stability is advisable. Unfortunately, precise 
mathematical theory does not exist to determine the 
optimal dissipation. Thus, most approaches reduce 
the approximation/flux/solution near the disconti- 
nuity to first order to achieve monotonicity and ro- 
bustness. Reduction to first-order accuracy locally 
results in the undesirable second-order 53 global ac- 
curacy if a uniform h - refinement is then performed. 

Artificial dissipation /filtering approaches are 
quite efficient and simple to implement in FD formu- 
lations. Unfortunately, they are often problem de- 
pendent, vary considerably on shock strengths, and 
are often user dependent. Although TVD and flux 
limiting approaches maintain monotonicity near dis- 
continuities, they unfortunately degenerate to first- 
order accuracy near smooth extrema. (See LeVeque 
88 for an overview of TVD techniques.) 

Essentially non oscillatory (ENO) 5 ‘ and later 
Weighted ENO (WENO) 90 - 74 schemes were de- 
veloped to circumvent nonlinear instability. ENO 
schemes choose the smoothest stencil from all avail- 
able design order stencils, thereby avoiding as much 
as possible interpolation/differentiation across dis- 
continuities. ENO schemes have been extremely suc- 
cessful algorithms over the past 10 years for prob- 
lems where both discontinuities and features requir- 
ing high-order spatial accuracy are required. See Shu 


129 for a detailed account of ENO/WEKO schemes 
and their applications. The principal difficulty with 
structured grid ENO schemes is their extension to 
complex geometries. The mathematical foundations 
for ENO/WENO schemes (bounded total variation 
proofs, etc.) are predominantly based on periodic 
or infinite domains, and are outside the context of 
boundaries. ENO/WENO schemes can not be im- 
plemented at several points next to boundaries be- 
cause they do not have smooth data outside the 
boundary to build high-order non-oscillatory sten- 
cils. The extension of ENO/WENO schemes to mul- 
tiple domains is complicated by numerous bound- 
ary interfaces throughout the domain. Another dif- 
ficulty with ENO schemes is that stencil searching 
algorithms have stencils that “switch” sometimes ar- 
bitrarily, which makes convergence to steady-state 
difficult. Atkins 5 proposed a smoothly varying 
stencil biasing technique to eliminate this problem. 
Integer stencil shifts were not allowed in the ap- 
proach. In addition, WENO schemes that are a 
smooth weighted sum of stencils in principle should 
not suffer from this difficulty. 

Durlofslcy, et. al. 4 ' and Abgrall 1 attempted to 
overcome the geometric complexity problems by 
building fully unstructured stencils using stencil 
searching algorithms. Ollivier-Gooch 108 suggested 
using a least-squares reconstruction approach for un- 
structured mesh ENO. Encouraging results were ob- 
tained for AGARD test case 1 with second- through 
fourth-order ENO, although convergence problems 
were experienced in the fourth-order case. 

The use of overset meshes is another approach 
used to overcome geometric difficulties. Wang and 
Huang 146 developed a compact ENO scheme and ap- 
plied it to a multi domain overset code. Boundary 
issues with the ENO formulation are still present, 
and additional complications of non conservation at 
the interfaces exist. Wanget al. 145 partially address 
the issue of interface conservation, but significant is- 
sues still remain for time-dependent discontinuous 
flows. 

Most high-order formulations are only guaranteed 
stable on linear problems, and some cannot even 
claim linear stability. The high-order DG FE formu- 
lations can claim a stronger form of stability. Barth 
14 and Barth and Chirrier 16 have designed numeri- 
cal fluxes that satisfy a nonlinear energy condition. 
They assume a convex entropy extension of the Eu- 
ler equations and bound the nonlinear “energy” of 
the system for all time in terms of the initial data. 
Simplified interface flux functions are derived to al- 
low this result. These results and others 66 provide 
an encouraging step towards nonlinearly stable for- 
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mulations that maintain high resolution. 

Although theoretically advantageous, high-order 
spatial discretizations in their present form still have 
several obstacles to overcome. Considerable work 
has been devoted to these methods over the past 
decade, yet few studies demonstrate the increased 
efficiency of high-order spatial discretizations for 
general-geometry aerodynamic simulations, includ- 
ing turbulence models. De Rango and Zingg 43 44 
address this specific question and achieve encourag- 
ing results for which high-order methods give more 
accurate solutions on a given grid. Additional work 
needs to be done to demonstrate increased efficiency 
for a given accuracy. 

Geuzaine et al. 52 and Delanaye et al. 41 apply 
the quadratic reconstruction FV scheme of Barth 
and Frederickson 12 to high Reynolds number flows. 
They achieve suitable convergence rates with gener- 
alized minimal residual (GMRES) and bi-conjugate 
gradient stabilized (Bi-CGSTAB) algorithms, and 
show second- and third-order convergence on irreg- 
ular and smooth meshes, respectively. They do not 
compare the efficiency of second- and third-order for- 
mulations, but note that the quadratic method con- 
verges more slowly than the comparable linear re- 
construction. Delanaye and Liu 42 report significant 
improvements in efficiency and accuracy in inviscid 
2-D calculation over a multi-element airfoil, compar- 
ing the quadratic and linear formulations. Results 
in three dimensions are not as dramatic. 

3.3.4 Langley Effort 

Langley has had a strong presence over the last 
decade in the following high-order spatial disciplines: 
1) structured grid FD, 2) structured grid ENO-FD 
and ENO-FV, 3) unstructured grid DG-FEM, and 
4) unstructured grid SUPG-FEM. The following in- 
formation is presented to summarize our experiences 
and provide guidance when comparing the different 
high-order methods. 

Table (1) can be used to compare the important 
attributes of current high-order formulations. Cat- 
egories are rated on a scale from one to five, with 
five being the best currently available, and one be- 
ing a capability representative of 1980. The cate- 
gories are 1) complex geometry, 2) grid adaptation, 
3) robustness (nonlinear), and 4) cost for a given 
accuracy requirement. At some level, all categories 
are closely related; however, we assume that each 
is independent from all others when assigning a nu- 
merical value. Specifically, the complex geometry 
category rates the capability of each method to ac- 
commodate complex 3-D configurations. This cat- 
egory is closely related to the locality of the dis- 
crete scheme. Grid adaptation is used to attack so- 


lution error. The second category, grid adaptation, 
describes each method’s success on adapted grids, 
including sensitivity to grid smoothness and ease of 
grid generation. The robustness category describes 
the robustness of the method for under-resolved fea- 
tures and discontinuities. In simple terms, this cat- 
egory rate whether the code “runs” (converges for 
steady-state cases, and does not diverge in time- 
dependent cases) with minimal user support. The 
cost category describes the cost of achieving a given 
accuracy. It is assumed that the necessary grid has 
been generated by whatever means are necessary in 
each case. 

The candidate schemes include two broad classes: 
1) unstructured database schemes, and 2) semi- 
st.ructured databases. We include the high-order 
FEM methods as unstructured methods because the 
structure within each element does not present a sig- 
nificant burden on the flexibility of the method. The 
unstructured candidate schemes are DG, SUPG, k- 
exact finite- volume, and k-ex act ENO FV. The can- 
didate structured grid schemes are Upwind FD, Up- 
wind FV, and WENO-FD. The numbers are sub- 
jective, and should only be used as a relative guide 
for the purpose of comparing strengths and weak- 
nesses. In general, researchers hold strongly varying 
opinions about the relative merits of each scheme. 

Table 1: Comparison of high-order schemes. Cat- 
egories are 1) complex geometry, 2) grid adaptation, 
3) nonlinear robustness, and 4) cost for a given ac- 
curacy. 


Method 

Geometry 

Adapt 

Robust 

Cost 






Unstructured 





DG FE 

5 

5 

3 

2 

SUPG FE 

5 

6 

2 

3 

A;- ex act FV 

5 

4 

2 

2 

LS-ENO FV 

4 

4 

4 

1 

Multi-block 





Upwind FD 

3 

2 

2 

5 

Upwind FV 

3 

2 

3 

3 

WENO-FD 

2 

1 

5 

3 


4 Convergence Acceleration 

A second algorithmic issue that contributes signif- 
icantly to the efficiency of the temporal algorithm is 
the convergence rate of the algebraic solution algo- 
rithm. At some point in the solution process the 
algebraic system 

A x = b (16) 
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is solved for the vector of data x, with 

I <9F 

A “ At + <9U 
x = U k+l -U k ; b = R(U k ) 


(17) 

(18) 


The matrix A is sparse. Equation (16) can be solved 
either directly, or iteratively. The dimensionality of 
x is 0(10') for 3-D problems, which makes direct 
methods uncompetitive with iterative solvers. Shur- 
complement methods can be used to attack equation 
(16) by subdividing A into numerous subproblems 
and interface conditions. Each subproblem is then 
solved directly, but the interface conditions are dif- 
ficult (dense). (See Saad 120 for details.) 

Iterative methods fall into two broad categories: 
stationary and nonstationary. 11 Stationary methods 
can be expressed in the simple form 


x k+1 = + c 


(19) 


with the matrix B independent of the iteration k. 
These methods are older, simple to implement, and 
usually not very effective when used alone. Sim- 
ple methods used in fluid mechanics are Jacobi, 
Gauss-Seidel (GS), symmetric GS, successive-over- 
relaxation (SOR), and symmetric SOR (SSOR). Nu- 
merous permutations of these methods exist includ- 
ing matrix reordering. (For details see Barrett et 
al. 11 ). More powerful stationary methods include in- 
complete factorization ILU(A), and block factoriza- 
tion. An ILU(A) approximately factors the original 
matrix A using an LU decomposition with the level 
of fill governed by the parameter Ic. Block factoriza- 
tion is motivated by tensor product grids (line data 
dependencies) where implicit solves along directions 
are efficient. Both are somewhat more expensive 
than basic stationary methods, but have consider- 
ably faster rates of convergence. Implementation of 
sweeping algorithms is complicated in a parallel en- 
vironment . 

Nonstationarv methods involve information that 
changes at every iteration. They are more recently 
developed, more difficult to understand, and more 
powerful. Commonly used examples in fluid me- 
chanics include conjugate gradient (CG), general- 
ized minimal residual (GMRES), bi-conjugate gra- 
dient (BiCG), quasi-minimal residual (QMR), and 
bi-conjugat-e gradient stabilized (Bi-CGStab). These 
methods update the solution in certain “directions” 
by considering inner products of current residuals 
and other Krylov space vectors arising during the 
course of the iteration. (See Saad 120 for details). 
On different problem classes, the convergence rate of 
nonstation ary methods varies considerably. Nacht.i- 
gal et al. 104 showed that a class of problems exists 


for which each of the aforementioned nonstationary 
methods is a clear winner (in terms of efficiency), 
and a clear loser (in terms of efficiency). 

A preconditioner is a matrix used to rotate a lin- 
ear system into a new system that has the same 
solution but is easier to solve in some sense. A 
left-preconditioner acts on equation (16) yielding the 
new system, 

M -1 A x = M -1 fc (20) 

where the new matrix M -1 A is easier to solve. 
The preconditioner changes the eigenstructure of the 
original system into a more compact set of eigenval- 
ues that an iterative method can attack more effec- 
tively. Ideally, a preconditioner should change the 
eigenstructure dramatically but at a minimal ad- 
ditional cost. Simple preconditioners used in fluid 
dynamics include the block Jacobi, GS, and SSOR 
methods. More powerful preconditioners include in- 
complete factorization ILU(A) and block factoriza- 
tion. 

Multigrid, at least in terms of elliptic problems, 
is a mechanism for rapid communication of mul- 
tiscale information. 23 Multigrid methods are usu- 
ally defined as a strategy to accelerate any sta- 
tionary or nonstationary iterative procedure. The 
solution is obtained on a sequence of grids, rang- 
ing from coarse to fine. Each grid smoothes the 
high-frequency components of the residual on that 
grid. Restrictions and prolongations communicate 
the data between grids. The resulting algorithm 
rapidly communicates long wavelength data via the 
grids, while damping short wavelength data by using 
efficient local smoothing operators. The exact choice 
of grid structure, restriction and prolongation oper- 
ators, and smoothers greatly influences the overall 
performance of the procedure. 

Multigrid methods are effective techniques for 
accelerating convergence of elliptic and hyperbolic 
problems. Convergence rates easily approach 0.1 
per cycle on elliptic problems such as the Poisson 
equation. The theoretical lower bound (if we rely' 
on coarse grid corrections) on the convergence for 
hyperbolic equations is 0.75 per cycle for a second- 
order spatial discretization, and is routinely achieved 
by general-purpose Euler solvers. (See Molder 103 
for details.) The convergence rate suffers consider- 
ably for high-Reynolds number (turbulent) viscous 
flow solutions. The primary cause for this slow- 
down is the highly stretched wall normal grids used 
to resolve turbulent boundary layers. Wall nor- 
mal stiffness can introduce stiffnesses of the order 
of 10 4 - 10 5 . A secondary cause of slowdown is the 
existence of significant regions of low Mach num- 
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ber flow around stagnation points and in recircula- 
tion regions. A third cause of slowdown comes from 
the imperfections of realistic grids. Generating grids 
that have 10 7 points without having grid anoma- 
lies is difficult. All general-purpose CFD solvers ad- 
dress each of these problems ultimately in the same 
way. For wall anisotropies, semi-coarsening and/or 
directional-implicit techniques precondition the wall 
normal boundary stiffness. Low' Mach number pre- 
conditioning is added in those regions of the flow 
below a critical Mach number. Grid anomalies are 
addressed with grid smoothing and movement algo- 
rithms in problem regions. 

Solution techniques within the aerodynamics com- 
munity are far from being “black-box” algorithms. 
Two decades of experience have shown that none 
of these algorithms is well suited for solving broad 
classes of high Reynolds number turbulent flows. 
Practitioners rely on combinations of a wide vari- 
ety of methods including 1) modified Newton-Krylov 
methods, 2) algebraic multigrid methods and 3) ge- 
ometric multigrid methods, 4) defect-correction it- 
eration techniques, and 5) sparse matrix methods. 

Anderson et al. 3 compared the efficiencies of 
several iterative strategies in the context of an 
unstructured, 3-D incompressible, Navier-Stokes 
solver. Multi element airfoils and high-lift sys- 
tems were used as test problems in 2-D and 3-D. 
The turbulence model used was that of Spalart and 
Allmaras. 132 GMRES w'as the Krylov method used 
in the study, with Gauss-Seidel or incomplete LU- 
decomposition used as a preconditioner. Newton- 
type solvers were shown to converge in the fewest it- 
erations. In terms of work and storage the multigrid 
algorithms are the most effective means of reducing 
the residual on the problems studied. 

In spite of all the powerful iterative techniques 
brought to bear on aerodynamic problems, the con- 
vergence rates for high- Reynolds problems can ap- 
proach 0.98 per cycle. Mavriplis demonstrated the 
capabilities of his unstructured code NSU3D 94 on 
complex 3-D configurations. 95, 96, 97, 98, 99, 100 The 
features in NSU3D are among the most advanced 
presently used in the CFD community. A wide va- 
riety of 3-D test problems was run including but 
not confined to 1) a realistic high-lift, configuration 
including a wing, pylon, and nacelle, 2) the trape- 
zoidal wing ", and 3) an ONERA M6. Grids ranged 
from 1-10 million vertices. Reynolds numbers were 
in the 1-10 million range with wall-normal spacing of 
10~ 5 - 10 -6 based on chord length. Published con- 
vergence rates for these cases ranged from 0.96 - 0.98 
per cycle. Approximately 500-1000 multigrid cycles 
were required to achieve residual levels converged to 


engineering tolerances. 

4.1 Bottlenecks 

The convergence rates quoted in the previous sub- 
section are based on steady-state solvers. The effi- 
ciency of time-dependent and steady computations 
are closely related, as the underlying nonlinear ma- 
trices are nearly identical. (The time-dependent for- 
mulation converges slightly faster for time-steps gov- 
erned solely by accuracy considerations). Thus, the 
cost of the unsteady calculations can be related to 
that of the steady. Simple back-of-the-envelope cal- 
culations reveal that at a minimum the unsteady 
computation will be equivalent to the solution of 100 
steady-state problems each having the same compu- 
tational complexity . Dropping an equation residual 
three orders of magnitude at a rate of 0.98 requires of 
approximately 400 iterations. Thus, modifying the 
convergence rate of the algebraic solution algorithm 
has a profound effect on the efficiency of the tem- 
poral algorithm. An order of magnitude increase in 
computational efficiency could be achieved by suc- 
cessful convergence acceleration efforts. 

The convergence characteristics of high-order spa- 
tial formulations has not been extensively studied. 
The dissipation level of the spatial operator fre- 
quently affects the convergence rate of the alge- 
braic system, with more dissipation producing faster 
convergence. High-order methods inherently have 
less dissipation, and could be more difficult to con- 
verge. An additional concern is the efficiency of 
multigrid methods as spatial order increases. The 
theoretical lower bound of multigrid methods (as- 
suming coarse grid corrections) on linear advection 
is ip > 1 — l/2 p where ip is the convergence rate, 
and p is the order of spatial approximation. Fourth- 
order methods should converge no faster than 0.9375 
per cycle. Ollivier-Gooch 108 used multigrid on 
an unstructured high-order FV scheme and expe- 
rienced increasing difficulty with convergence as the 
order of the approximation increased. Delanaye et 
al. 41 notes that quadratic elements converge more 
slowly than equivalent linear elements, when a differ- 
ence GMRES algorithm is used. Interestingly, Bon- 
haus 21 using SUPG with a Newton-Krylov GMRES 
method and diagonal preconditioning, experienced 
no changes in convergence rate with increasing or- 
der. 

4.2 Langley Effort 

Current iterative nonlinear solvers require 0( 10 3 ) 
iterations to converge, based on high Reynolds num- 
ber (ps 10 b — 10‘), turbulent, separated, 3-D, com- 
plex geometry flows. Newton solvers converge these 
problems in C4 ( 1 0 1 ) iterations if a reasonable ini- 
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tial guess is given. Currently, a. group at Lang- 
ley is studying methods that potentially have text- 
book multigrid efficiency (TME). 22, 23 Methods that 
achieve TME converge at rates that are independent 
of the number of degrees of freedom (grid), and con- 
verge to the truncation error of the discretization 
in approximately 10 1 work units. A work unit is 
defined as the work equivalent to the evaluation of 
the residual. Thus, TME methods have a poten- 
tial increase in efficiency of 0 ( 10 2 ) compared with 
existing state of the art solvers. The potential effi- 
ciency of TME methods relies the “factorizability” 
property of the Navier- Stokes (NS) equations. Split- 
ting the NS equations into factors decouples the sys- 
tem somewhat, and allows optimal (fast and inex- 
pensive) operators to be constructed for each por- 
tion. This “divide and conquer” approach yields 
nearly optimal efficiency for the entire problem. The 
Langley effort was showcased at the 2001 AIAA 
CFD conference. 24 ! 4Si 119, 13 °- 139 Work continues on 
this revolutionary method to implement TME algo- 
rithms on general purpose aerodynamic solvers. 

5 Turbulence Modeling 

5.1 Overview: Current practices 

An equally important aspect of temporal algo- 
rithms is the underlying turbulence models being 
solved. Numerous turbulence models exist for at- 
tached steady turbulent flows. During the last two 
decades great strides have been made in tuning these 
turbulence models to increase their robustness and 
generality. Practitioners routinely utilize these mod- 
els to predict the flow behavior of a surprisingly 
broad class of complex problems with acceptable 
confidence levels in their solutions. Unfortunately, 
the same level of maturity does not exist for nonsta- 
tionary flows, where time-averaged turbulence quan- 
tities inadequately describe the important dynamics 
of the flow. Oftentimes, these flows are dominated 
by large-scale flow features that are not properly 
modeled by conventional turbulence models. Flows 
with massive separations such as bluff-bodv wakes, 
cavity flows, shock-induced separations, and recircu- 
lation zones almost always fall into this category. 

5.1.1 LES 

In 1970 Deardorff 39 published the first results of 
a large-eddy simulation (LES). The objective of an 
LES is to simulate, or directly compute, the large 
energy-containing fluid motions and to model only 
the small scales that are unresolved by the grid (the 
subgrid scales). A subgrid-scale (SGS) model acts to 
remove the energy associated with the small scales 
and therefore facilitates the global energy transfer 


from the large scales to the small scales. When prop- 
erly implemented, LES can be used to simulate the 
turbulence in a flow at low to moderate Reynolds 
numbers. For more complete discussions of the con- 
cepts and applications of LES see Pionrelli 115 ' 116 
and the references therein. 

In an LES calculation, the smallest resolved scales 
are determined by the grid-cell size. On finer grids, 
more of the flow is simulated and less is modeled. 
The SGS model therefore has an explicit dependence 
on the local grid-cell size. This feature of LES has 
led some to believe that the grid-cell size is unre- 
stricted and simply corresponds to the break in re- 
solved and subgrid scales. This sort of thinking usu- 
ally leads to poor calculations. An important as- 
sumption in aii LES calculation is that the energy- 
containing scales are actually simulated. To do this, 
the peak in the energy spectrum must be in the re- 
solved range of scales and the cutoff between the 
resolved and subgrid scales should be in the inertial 
wavenumber range. 

The development of dynamic subgrid-scale mod- 
els for LES has motivated the use of LES on non- 
stationary turbulent flows over extremely complex 
configurations. However, estimates of the grid re- 
quirements for LES computations over realistic ge- 
ometries at realistic Reynolds numbers indicate that 
LES is not likely to be a viable option for most flows 
for several decades to come. 134 

5.1.2 URANS 

Despite the abovementioned shortcomings of ex- 
isting turbulence models, most time-dependent 
codes use unsteady extensions of popular steady- 
state algorithms. These calculations are typically 
referred to as unsteady Reynolds-averaged Navier 
Stokes (URANS). The turbulence models used in 
URANS span the spectrum of available models. 

Hold et al. 65 used a Baldwin-Lomax model to 
compute the unsteady flow about a rounded half- 
cylinder embedded in a flat plate. Their time- 
averaged pressure coefficients agreed reasonably well 
with experimental measurements. Surface pressure 
spectra, however, showed far more disagreement 
with the measured spectra. Hold et al. 65 suggested 
that the disagreements were primarily temporal and 
spatial resolution issues rather than a turbulence 
modeling problem. Their paper illustrates a signif- 
icant problem associated with assessing turbulence 
models for unsteady flows; unless care is taken to 
first, assess grid and time-step issues, the influence of 
the turbulence model is almost impossible to evalu- 
ate. Because of the computational expense of doing 
thorough resolution checks for unsteady problems, 
these numerical issues are often left unresolved. 
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Khorrami, Berkman, and Choudhari 79 and later 
Khorrami, Singer, and Berkman 80 used the two- 
equation SST (k-u;) model of Menter 101 for 2-D 
unsteady calculations of flow in the vicinity of a. 
leading-edge slat. In their work the grid resolution 
and time-step were sufficiently fine for them to sus- 
pect that the turbulence model was responsible for 
excessive damping of the large eddies in the slat 
cove. Subsequent work by Khorrami, Singer, and 
Lockard 81 locally eliminated the turbulence produc- 
tion term in much of the slat cove and observed un- 
steady coherent vortex motions that closely resem- 
bled the PIV measurements reported in the work of 
others. 112 - 140 - 141 

Scotti and Piomelli 121 assessed four low Reynolds 
number turbulence models in pulsating channel flows 
generated by oscillating the imposed pressure gra- 
dient. The apparent simplicity of these flows was 
deceptive. In fact the pulsating flows represented 
significant challenges to the URANS models. The 
turbulence was out of equilibrium, in some cases 
so much so that partial relaminarization occurred, 
followed by a re-transition. Somewhat surpris- 
ingly, the phase- averaged streamwise velocity pro- 
files matched corresponding LES results reasonably 
well when plotted in wall coordinates. However 
the time-averaged velocities of the URANS mod- 
els were always less than those determined from the 
LES. Significant differences in the Reynolds shear 
stresses and the turbulent kinetic energies were also 
observed. 

5.1.3 Composite LES/RANS Schemes 

Recently, a class of composite models has been de- 
veloped for unsteady flows. These composite models 
attempt to blend the unsteady capabilities of LES 
with a method having grid requirements that are 
more like those conventionally used in Reynolds- 
averaged Navier-Stokes (RANS) calculations. A 
composite model will typically involve a RANS tur- 
bulence model and RANS- type grid in regions near 
solid surfaces, where the resolution of turbulent ed- 
dies would require exceptionally fine grid resolution. 
In these regions typical RANS models are usually 
adequate for modeling the effects of turbulent flow. 
Away from wall regions, where large unsteady flow 
structures dominate the flow, the composite model 
shifts from a RANS- type turbulence model to a less- 
dissipative LES- type subgrid-scale model. Similarly, 
away from wall regions, the gridding strategy will 
shift from a RANS- type grid to a grid amenable to 
resolving at least some of the unsteady turbulence 
scales of motions. 

A variety of different composite models have been 
proposed in the last few years. The composite 


models go by different names and involve different 
turbulence models for the RANS and LES regions 
and different switching functions to shift between 
the RANS and LES regions. Batten, Goldberg, 
and Chakravarthy 18 proposed limited Navier Stokes; 
Zhang, Bachman, and Fasel 148 developed their flow- 
simulation methodology; Arunajatesan, Shipman, 
and Sinha 4 developed the hybrid RANS-LES. How- 
ever, the first, and perhaps the most thoroughly 
tested of these composite models has been the de- 
tached eddy simulation (DES) model of Spalart. 135 
In its most common form, the DES model is imple- 
mented in combination with the Spalart- Allmaras 
(SA) turbulence model, 133 although recently DES 
results with Menter ’s SST model 101 have also been 
published. 138 In either case, the switch from the 
RANS region to the LES region is governed by a ra- 
tio of grid-cell size with distance from the wall. Vis- 
cous grid cells (high-aspect ratio with lengths long 
compared to the distance from the wall) characterize 
the RANS regions. Small, approximately uniformly 
sized grid cells far from solid surfaces characterize 
the LES region. In theory, the RANS portion cor- 
rectly models the rather benign attached turbulent 
flow regions, while the LES portion has sufficient 
grid resolution and sufficiently reduced eddy viscosi- 
ties to simulate the large-scale turbulent flow struc- 
tures. 

The ease with which DES can be implemented in 
a code is a mixed blessing. Upgrading an existing 
SA model to accommodate the DES model requires 
minimal additional logic. The triviality of the modi- 
fications required for this transformation has led nu- 
merous well-intentioned users to perform and pub- 
lish DES-like calculations. Unfortunately, it is easier 
to obtain results than to verify and critically inter- 
pret them. 

Proper grid resolution is always important in 
CFD. However, with DES, the grid takes on new- 
significance, as the grid determines the sw'itch be- 
tween the RANS region and the LES region. For- 
tunately, Spalart 136 provides some guidance for the 
development of appropriate grids. Importantly, once 
the RANS region of the flow is adequately resolved, 
improved grid resolution should be undertaken lo- 
cal to the LES portion of the simulation, not glob- 
ally. Over-resolving the RANS regions can displace 
the RANS/LES switch-over surface into the bound- 
ary layer, where the grid is likely to be insufficient 
for a reliable LES calculation. Of course, continued 
grid refinement will eventually lead to adequate LES 
grids over the entire flow, although for aerodynamic 
flows of commercial interest, solving the flow- on such 
a grid will likely be a decades-long endeavor. 
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With a DES grid well designed for a particular ap- 
plication, global grid refinement is likely to adversely 
affects the results. Nikitin et al 105 discusses these 
circumstances in the context of a series of channel- 
flow calculations. Although modifying the switching 
function, as is done in other composite LES/RANS 
methods, can alleviate some of the deleterious ef- 
fects, a key to successful implementation of these 
methods is going to involve careful attention to grid- 
ding. 

Another concern with the implementation of com- 
posite methods is the numerical diffusion, especially 
in the LES region where diffusion associated with the 
model is far less than in the RANS region. In simula- 
tions of the unsteady flow about a sphere, Constan- 
tinescu and Squires 38 used both second-order and 
fifth-order accurate upwind differences for the con- 
vective terms. All other operators were discretized 
using three-point, second-order accurate central dif- 
ferences. On the same grid, they found that the 
results of calculations with the fifth-order accurate 
scheme agreed better with a separate LES calcula- 
tion than the results obtained with the second-order 
upwind scheme. 

In later work, Strelets 138 was concerned about the 
diffusion of upwind schemes in the LES region and 
the lack of stability of central-difference schemes in 
the RANS regions. To weaken the adverse effects 
of upwinding in the LES regions, he used a hybrid 
central/upwind approximation of the inviscid fluxes. 
The hybrid is designed so that in the RANS region 
the scheme is “almost upwind” and in the LES re- 
gions the scheme is “almost centered.” 

Recent work at NASA Langley suggests that 
second-order central-difference schemes can be used 
for both the RANS and LES regions, although the 
grid resolution must be chosen appropriately. The 
key is to ensure that the turbulence model, and not 
the numerics, controls the diffusion. 

DES has been used in unstructured grid environ- 
ments, first by Forsythe, Hoffmann, and Dietiker. 
Their research suggested that Coes . the one addi- 
tional free parameter in DES (as compared to SA) 
needed to be adjusted when DES is used in conjunc- 
tion with tetrahedrons. Because Coes weights a 
measure of grid-cell size in the RANS/LES switch- 
ing function, the need to adjust Coes when differ- 
ently shaped grid elements are used is not surprising. 
Pelaez, Mavriplis, and Kandil 113 used the standard 
value of Coes for their unstructured grid calcula- 
tions and ultimately concluded that the value should 
be determined by performing a decaying homoge- 
neous turbulence test case. 


In spite of the many unanswered questions associ- 
ated with DES, in the right hands it can be a useful 
tool for simulating unsteady complex flows. Travin 
et al. 143 used DES for flows around a circular cylin- 
der. They explored the effects of grid, Reynolds 
number, 2-D versus 3-D. DES versus LIRANS, and 
laminar versus turbulent separation. An important 
conclusion, first suggested by Spalart et al. 134 is that 
2-D DES calculations are not fruitful. To achieve the 
advantages of DES, a true LES should be done in 
the LES regions which implies that the full 3-D flow' 
needs to be simulated. Without the third dimension, 
important turbulence dynamics are inoperative. An- 
other important conclusion is that the benefits of 
DES vary with the flow and the information desired. 
In the case of a circular cylinder with laminar sepa- 
ration, DES proved itself superior to URANS, even 
for obtaining time-averaged quantities. However, for 
the case of turbulent separation, the time-averaged 
quantities were the same, regardless of whether DES 
or URANS was used. Hedges, Travin, and Spalart 39 
found that URANS performed almost as well as DES 
for time-averaged quantities on a four-wheel landing- 
gear model. The unsteady flow around the landing 
gear appeared much more turbulent-like in the case 
of DES compared to URANS, but unfortunately un- 
steady flow data was not available for comparison. 

5.2 Bottlenecks 

In this section we will review some of the bottle- 
necks that are holding back further progress in re- 
solving turbulence modeling issues. In particular, we 
will discuss experimental validation problems, con- 
sistency concerns for the composite methods, and 
the high cost of the calculations. 

5.2.1 Validation by Experiments 

One of the great difficulties with new turbulence 
models for unsteady flows is the lack of relevant ex- 
perimental data with which to validate the mod- 
els. As discussed earlier, different types of turbu- 
lence models with very different unsteady character- 
istics can provide very similar time-averaged results. 
Hence, if the time-dependent behavior of the flow is 
important in an application, then experiments that 
provide time-dependent data are indispensable. 

Although flow past circular cylinders and spheres 
has provided useful validation tests in the past, for 
applications that involve flow actuators, particu- 
larly those that involve unsteady mass addition, new 
types of experiments are required. 

A typical unsteady flow actuator comprises a noz- 
zle or orifice that communicates with a cavity. The 
cavity includes either a flexible membrane or a mov- 
ing wall. Unsteady deformation of the flexible mem- 
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ber inside the cavity results in a local pressure 
change that drives fluid through the nozzle and into 
an external flow. In practical applications, the ex- 
ternal flow is turbulent without the presence of the 
actuator; the flow inside the actuator is often a mix 
of laminar and turbulent regions. 

A good place to start with the experiments is 
where a number of calculations of flow actuators 
have begun, i.e. , with actuators in no-flow and lam- 
inar flow environments. Such experiments are cur- 
rently being performed at Langley. These experi- 
ments are time-consuming and require attention to 
details not typically recognized as important. For in- 
stance, the deformation of the flexible member may 
need to be simulated or modeled with greater fidelity 
than expected, time lags between the electrical ac- 
tuation of the flexible member and its actual me- 
chanical actuation need to addressed, and flow mea- 
surements in largely inaccessible locations, like the 
actuator cavity, are desired. 

The problems prove much more difficult when the 
external flow is turbulent. In that case, character- 
istics of the external turbulent boundary layer must 
be measured. In addition, the interaction of the flow 
structures in the turbulent boundary layer with the 
flow in the cavity must be characterized. Validated 
criteria for choosing a particular turbulence model- 
ing approach for a specific type of actuator flow have 
not been established. However, we speculate that if 
a significant amount of the external flow is ingested 
into the actuator cavity, then a full LES calculation 
may be the only reasonable approach. On the other 
hand, if the actuator is biased such that it does not 
ingest the external flow, a DES calculation might be 
able to calculate the relevant flow features. One con- 
cern with DES in this sort of application is that the 
relevant unsteady flow regions are close to the wall; 
hence the model may revert to its RANS character- 
istics in these regions and be too diffusive. Other 
composite LES/RANS models may be more appro- 
priate for such flows. Although explorations of these 
problems are ongoing at Langley, no definitive guid- 
ance can be provided to the CFD community yet. 

5.2.2 Consistency Concerns 

Gatski 50 raises questions as to the formal con- 
sistency of the composite methods. The compos- 
ite methods assume that the flow variables from the 
RANS regions match smoothly with the correspond- 
ing variables from the LES regions. Gatski 50 shows 
that such a match cannot be taken for granted, even 
without considering the consistency issues associ- 
ated with the coupling of the temporal averages used 
to derive RANS models and the spatial filtering used 
to derive SGS models for the LES regions. To date, 


these consistency issues have not been addressed by 
the modeling community. 

5.2.3 Calculation Cost 

One of the problems with using a composite 
LES/RANS method is that for the LES region to 
be reasonably simulated, it needs to be not only 
time-dependent, but also three-dimensional. With- 
out three-dimensionality, the vortex stretching and 
tilting mechanisms, which are so important in real 
flows, cannot be appropriately simulated. The ne- 
cessity to do 3-D calculations dramatically raises the 
costs of the calculations, both in terms of computer 
memory and in terms of run time. Turnaround time 
for circular cylinder calculations on eight processors 
of an SGI Origin 2000 can easily be multiple weeks. 
Grid and time-step studies therefore turn the activ- 
ity into a very long process, even for such a relatively 
simple flow. The time can be significantly reduced 
through the use of more extensive parallel process- 
ing, but the parallel efficiency must be maintained 
over approximately 100 processors for the calcula- 
tion turnaround times to become sufficiently short 
for good systematic studies to be performed. Those 
like Travin et al. 143 who have performed such studies 
have added greatly to our understanding of the per- 
formance of these methods. More work along those 
lines needs to follow. 

5.3 Langley Effort 

The efforts at Langley in this area have concen- 
trated on the development of appropriate validation 
experiments and also in doing some of the careful 
studies needed to calibrate the performance of DES 
with different codes. Some additional developmen- 
tal work on new families of composite LES/RANS 
models has been funded by NASA Langley Research 
Center. Although the DES model has by far the 
greatest, following, and therefore the most thorough 
testing, some aspects of the model raise questions 
as t.o its suitability for actuator-type computations. 
In particular, the fact that the boundary between 
the RANS region and the LES region is entirely dic- 
tated by grid size and distance from the wall makes 
the simulation of near-wall recirculation regions dif- 
ficult. Other versions of composite models may work 
better in flows where near-wall recirculation regions 
are anticipated. 

Another important role played by Langley is to 
act as a clearinghouse for ideas related to unsteady 
turbulence models. Langley researchers can provide 
appropriate suggestions and criticism of emerging 
concepts so that the developers can address these 
concerns early in the models’ development. As an 
example of this interaction, the consistency concerns 
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discussed earlier are now being addressed in at least 
one new model under development. 

At NASA Langley Research Center a concerted 
effort is underway to measure and compute the flow 
in active flow control devices. This effort brings to- 
gether experimentalists, turbulence modelers, and 
computationalists who all have been working to- 
gether to develop appropriate techniques for the sim- 
ulation of actuator flows. 

6 Actuator Boundary Conditions 

6.1 Overview: Current Practices 

Interest in active flow control for drag or noise 
reduction, flow vectoring, mixing enhancement, and 
separation control has stimulated the recent develop- 
ment of innovative synthetic jet actuators that cre- 
ate localized disturbances in a flowfield. Synthetic 
jets are generated by a dynamic fluid actuator con- 
sisting of a cavity enclosed by one (or more) moving 
diaphragm (s) driven into transverse oscillations at 
their resonance frequency. The distinctive feature 
of these actuators is that they have minimal power 
requirements and have jet-like characteristics with- 
out the need for mass injection. Although net mass 
injection into the overall system occurs over each 
cycle of operation, the momentum transfer into the 
embedding flow is nonzero. These features enable 
synthetic jets to effect significant global modifica- 
tions in the embedding flow on scales that are one 
to two orders of magnitude larger than the charac- 
teristic length or force scale of the jets themselves. 

In recent 

years, because of considerable improvements in com- 
puter resources, more attention is being devoted to 
numerical simulation and optimization of synthetic 
jet actuators. 82 ' 46 58 - 33 - 147, 29 - 64 - 89 - 118 7 " 87 The 
presence of several temporal and spatial scales and 
moving boundaries in the problem makes simula- 
tion of such unsteady flows computationally very ex- 
pensive. To reduce the computational cost, several 
approaches have been developed. All the methods 
can be divided into two classes: 1) simplified models 
without simulation of the flow' inside the actuator, 
and 2) numerical simulation of both the exterior and 
cavity flows. 

In the first class of methods, a synthetic jet gener- 
ated by harmonic motion of the actuator diaphragm 
is modeled by using simplified boundary conditions 
imposed at the orifice exit. In the work of Donovan, 
Krai, and Cary, 82 and Krai et al. 4S , the flow within 
the cavity is not calculated and the perturbation to 
the flowfield is introduced through the wall-normal 


component of velocity at the orifice exit 

v{x,y = 0 ,t) = VA(x)sin(ujt) 
u(x, y = 0,f ) = 0 


where x and y are the flow streamwise and cross- 
stream directions, respectively. Different spatial 
variations of Va(x) over the orifice have been con- 
sidered. Numerical experiments 4S> 82 indicate that 
a “top hat” distribution most closely matches the 
experimental data. A modified boundary condition 
on the pressure at the orifice is introduced through 
a consideration of the normal momentum equation. 
Taking into account time-harmonic velocity pertur- 
bations, this condition, obtained under the assump- 
tion that the flow is incompressible, becomes 

dv 

— = —pVA(x)uicos(ut) (22) 

dy 

The numerical simulations based on the boundary 
condition (eqs. 21 and 22) show good qualita- 
tive agreement with the experiments. 123 ' 122 ’ 131 A 
similar approach modeling the actuator as a blow- 
ing/suction boundary condition, which can be fully 
specified in advance of the calculation, is used by 
others. 29, 5S - ® 4, 89 

An alternative technique for modeling synthetic 
jet actuators has been proposed by Carpenter et al. 33 
This theoretical model for the actuator is based on 
classic thin plate theory for the diaphragm dynam- 
ics. The flow through the orifice is modeled using 
unsteady pipe-flow theory. This approach is based 
on the assumption that streamlines in the orifice exit 
are parallel to its axis, which is an adequate approxi- 
mation if the length-to-diameter ratio is much larger 
than unity. The governing equation for the orifice 
flow is given by 

dv pv\v\ dp d ( 1 dv\ 

P ~di + ~W = + [rfr) 

where y and r are the axial and radial coordinates, 
respectively, v is the axial velocity, l is the orifice 
length, and density p is taken to be the instanta- 
neous mean of the cavity and external densities. The 
inertial term pvv y is modeled approximately by the 
second term on the left-hand side of equation (23). 
The dynamics of the air in the cavity is ignored and 
the pressure there is calculated by means of the per- 
fect gas law. The cavity and the external boundary 
layer flowfields are linked by requiring continuity of 
velocity and pressure at the orifice exit. In Carpen- 
ter, Lockerby and Davies, 33 only the blowing phase 
of the actuator dynamics has been studied. 
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The second class of methods is based on a. direct 
numerical simulation of the entire problem includ- 
ing the flow inside the actuator. Rizzetta, et al. 118 
solve the unsteady compressible Navier-Stokes equa- 
tions in the external region, the cavity itself, and 
the throat of the actuator on separate grids that are 
linked with each other through a Chimera method- 
ology. The membrane motion is simulated by vary- 
ing the position of appropriate boundary points. As 
follows from the numerical calculations presented in 
reference [ 118 ], the internal cavity flow becomes peri- 
odic after several cycles. Therefore the velocity pro- 
file across the jet exit at each time-step was recorded 
and was used as a boundary condition in subsequent 
runs involving the external domain only. For com- 
putations that consider only the upper exterior do- 
main, the transverse and span- wise velocity compo- 
nents (orthogonal to jet axis) are set to zero and the 
inviscid normal momentum equation, expressed as 


dp ( dv dv\ 

dy P \dt + V dy) 


(24) 


is used to establish the pressure. The orifice exit 
density is extrapolated from the interior solution. 
This approach provides more accurate description 
of the flow details at the orifice than the simplified 
boundary condition of equations (21) and (22). Sim- 
ilar direct numerical simulation of the external and 
cavity flows has been performed by Joslin et al. 77 
and shown good agreement with experimental re- 
sults. 

To avoid the integration of the Navier-Stokes 
equations on a moving grid, an alternative technique 
is used by Lee and Goldstein. 87 The main idea of the 
method is to impose a localized body force along de- 
sired points in the computational mesh to bring the 
fluid there to a specified velocity so that the force has 
the same effect as a solid boundary. The desired ve- 
locity is incorporated in an iterative feedback loop 
to determine the appropriate force. For a moving 
boundary with velocity Vj^aqt), an expression for 
the body force is 


t 

F(x, t) = a j (v - VA)dt' + fi(v - Va) (25) 
o 

where v is the fluid velocity, and a and /? are 
user-defined constants that are negative and can be 
treated as the gain and damping of the force field 
with dimensions of M / (L 3 T 2 ) and M/L 3 T, respec- 
tively. This approach allows a moving diaphragm 
without using a time-dependent coordinate trans- 
formation. 


6.2 Bottlenecks 

In spite of the fact that the methods mentioned 
above have successfully been used for modeling syn- 
thetic jet actuators, several issues persist. One of 
the main drawbacks of the first class of methods is 
that the simplified boundary conditions do not pro- 
vide conservation of mass, momentum, and energy 
through the actuator orifice. Because these methods 
use the normal momentum equation to calculate the 
pressure, whereas the other quantities are extrapo- 
lated or prescribed analytically, this boundary con- 
dition does not satisfy the governing equations at the 
boundary and, therefore, does not provide the con- 
servation properties. Another disadvantage of the 
boundary condition [eqs. (21) and (22)] is its in- 
ability to account for changes in the pressure field 
caused by the interaction of the external boundary 
layer and actuator. Furthermore, as has been shown 
in Lee and Goldstein, 87 the real streamwise velocity 
profile and the velocity component in the cross-flow' 
direction are far from the analytical expressions of 
equation (21). 

The main problem associated with the second 
class of methods is complexity /cost. The numerical 
calculation of the cavity flow requires a large number 
of grid points. For geometries fitted w'ith multiple 
actuators, grid requirements for the actuators could 
exceed those of the exterior flow, and would con- 
tribute extensively to the computational cost. An- 
other consideration is the actuator Mach number, 
which varies from 0.001 (near the diaphragm) to 0.1 
(at the orifice exit). This variation of the flow pa- 
rameters from fully incompressible in the actuator 
to fully compressible in the exterior region imposes 
very severe requirements on a numerical method and 
increases the algorithm complexity. 

6.3 Langley Effort 

As follows from the literature overview presented 
above, the research in the area of active flow control 
is of empirical nature, mainly due to the computa- 
tional cost involved and lack of confidence in compu- 
tational methods for such complex time-dependent 
flows. A strong effort 141 is currently underway 
at. Langley toward constructing a new methodology 
that combines the accuracy and conservation prop- 
erties of the simulation methods with the efficiency 
of the techniques based on simplified boundary con- 
ditions. In contrast to the methods found in the 
literature, the new approach uses a reduced-order 
model of 2-D or 3-D actuators. In other words, 
the multidimensional actuator is simulated by solv- 
ing the time-dependent 1-D Euler equations similar 
to those used for the quasi-one-dimensional nozzle 
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environment . 


blowing/suction boundary condition 



Time 


Figure 3. Comparison of the mass rate errors 
for the 2-D actuator problem obtained with the 
quasi-l-D actuator model and the nonconservative 
boundary condition (21,22). 


problem. The simplified actuator model has several 
advantages. First, this approach provides conser- 
vation of not only mass, but also momentum and 
energy. Second, the new method is much more effi- 
cient in terms of computational time compared with 
the 2-D or 3-D numerical simulation of the flowfield 
in the actuator. Third, this reduced-order model 
retains some multidimensional features of the realis- 
tic actuator that are governed by the area variation 
of the quasi-one-dimensional nozzle, and, therefore, 
can be used for qualitative study of the resonance 
characteristics of the actuator. As follows from the 
comparison presented in figure (3). the conventional 
boundary conditions based on the normal momen- 
tum equation do not provide mass conservation. The 
maximum mass rate error, which occurs during the 
suction stage, is of the order of 15 percent if the 
normal momentum equation is used as a boundary 
condition for pressure. In contrast to the conven- 
tional approach, the new method provides conserva- 
tion of mass, momentum, and energy. As a result, 
the mass rate error is reduced by one to two orders 
of magnitude compared with that obtained with the 
blowing/suction boundary condition [eqs. (21) and 
(22)]. These preliminary results are very encourag- 
ing. In our future work, we will focus on calibration 
of the new methodology with readily available ex- 
perimental data and numerical simulation of time- 
dependent. flows encountered in active flow control 


7 Conclusions 

The current status of time-dependent algorithms 
is presented. Special attention is given to algo- 
rithms used to predict fluid actuator flowfields. The 
overview begins by considering algorithmic issues 
that could greatly improve the temporal efficiency 
of actuator simulations. The general state of time- 
dependent. turbulent models for nonstationary flows 
is then assessed. Finally, an efficient, new fluid actu- 
ator boundary condition is presented. Each section 
begins by describing the current state of the art in- 
cluding notable impediments in the field, and con- 
cludes with a summary of current. Langley efforts. 

Profound improvements in the efficiency of tempo- 
ral algorithms could be achieved in the next, decade. 
Notable leverage in time-dependent, methods exists 
in the following algorithmic areas: 1) implementa- 
tion of high-order (p > 3) temporal schemes, 2) 
implementation of high-order ( p > 3) spatial algo- 
rithms, and 3) convergence acceleration techniques 
for complex high-Reynolds number flows. Signifi- 
cant, impediments exist in each of these three cate- 
gories. 

High -order schemes need huge time steps to uti- 
lize their full potential. Algebraic solvers with rapid, 
time-step independent, convergence rates are neces- 
sary for these schemes. The principal impediment 
facing the implementation of high-order temporal 
schemes is the need for robust and rapid algebraic 
equation solvers. Error estimation, variable time- 
stepping and automated iteration termination will 
immediately follow. Current fourth-order schemes 
are asymptotically close to being optimal, and fur- 
ther increases in efficiency are difficult to obtain. 

High-order spatial operators must be flexible 
enough to accommodate complex geometries, grid 
adaptation and nonlinear instability. Methods uti- 
lizing compact locality (unstructured methods) are 
advantageous when addressing complex geometries 
and grid adaptation. Unstructured high-order for- 
mulations include Finite-Elements (FE), fc-exact Up- 
wind Finite- Volume (FV), or fc-exact ENO FV. The 
principal impediment facing the implementation of 
high-order spatial operators is their lack of nonlinear 
stability resulting from unresolved features and dis- 
continuities. Extensive research is currently under- 
way within the FE community. Some formulations 
possess a nonlinear L 2 stability property. It is still 
unknown whether this stability is a strong enough 
for general purpose solvers, or whether stronger sta- 
bility conditions (TVD, TVB, ENO) must be pur- 
sued. 
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A wide variety of different iterative methods are 
currently used. Nevertheless, the convergence rate 
of current state of the art iterative solvers is poor. 
Complex high Reynolds number 3-D simulations 
converge at per cycle rates in the range 0.96 — 0.98. 
The principal impediments for rapid convergence are 
high aspect ratio cells in the turbulent boundary lay- 
ers, and low Mach number regions in the flow. 

Conventional URANS turbulence models are not 
very accurate for nonstationary turbulent flows. A 
new class of composite LES/RANS schemes has been 
developed to address these inadequacies. Among 
these is the detach eddy simulation (DES) approach, 
and is currently being verified and validated. The 
DES scheme is very expensive because is requires 
a full 3-D time-dependent simulation of the flow in 
question. The principal impediment facing all com- 
posite approaches is validation against high quality 
experimental data. The accuracy and consistency of 
the blended region between the outer LES and inner 
RANS regions is a concern with composite methods. 
The accuracy of this region is extremely important 
in fluid actuator simulations. 

An optimistic estimate indicates that successful 
research could improve overall efficiency by O(10 1/,J ) 
for temporal algorithms, by O(10 3 ' 2 ) for spatial al- 
gorithms, and by O(10 3 ' 2 ) for convergence accelera- 
tion. Improvements in the modeling of fluid actua- 
tors could account for a factor of two. Moore's law 
predicts that increases in computer hardware will 
yield O(K) 1 ) improvements. All these effects can be 
combined in a multiplicative sense to yield poten- 
tial improvements of O(10 4 ). Improvements of this 
magnitude would allow us to do 3-D optimization 
studies based on DES turbulence models, including 
fluid actuator design and resonance, and actuator 
placement and coupling studies. 
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